clc;
clear all;
close all;

tic % STOPWATCH ON

inter_len = 600;       %交织长度，300ms或者600ms
inter_tp = 2;       % 0：不交织； 1：行列交织； 2：素数交织
chann_flag = 1;
fig_on = 0;
no_noise = 0;
iter_max = 4;
if(inter_len==600)
    inter_p = 133;           %素数交织的素数
    inter_line = 72;        %行列交织的行数
else
    inter_p = 67;           %素数交织的素数
    inter_line = 36;        %行列交织的行数
end
fid = fopen('result.txt', 'a');
str=datestr(clock);
fprintf(fid,'\n\n*********************%s*********************start\n', str);
% =========================================================================
sps = 8;
Rb = 4000;
Tb = 1 / Rb; % bit duration
Ts = Tb*(1/sps); 
Fs = 1 / Ts;
BW = 6000;
BW_bb = BW / 2;         % 基带带宽

%% Pulse shape & Variable initialization
pulse            = 2;      % 1 -> lorentzian pulse
                           % 2 -> GMSK pulse BT = 0.3
                           % 3 -> LRC pulse
                           % 4 -> LREC pulse
pulse_length     = 4;      % 1  -> Full response
                           % >1 -> Partial response
M_ary            = 2^1;    % M_ary symbols used 2 -> Binary
h_mod            = 0.5;    % Modulation index
[g_t,q_t]     = pulse_gen(pulse, pulse_length, 1, Fs, sps);


d_lpf = fdesign.lowpass('Fp,Fst,Ap,Ast', 0.375, 0.45, 0.1, 60);
% d_lpf = fdesign.lowpass('N,Fp,Fst',51, 0.1875, 0.3);
Hd_lpf = design(d_lpf);
% fvtool(Hd_lpf);
h_lpf = Hd_lpf.Numerator;
odr_lpf = (length(h_lpf)-1) / 2;

% load 'sh_1d3_24_128.mat';
% load 'sh_1d5_16X20_msk.mat';
if (inter_len==600)
    load 'sh_1d3_12X18_60_msk.mat';     %按照0.6s交织，则这改成30
    z = 60;         %按照0.6s交织，则这改成30
else
    load 'sh_1d3_12X18_30_msk.mat';     %按照0.6s交织，则这改成30
    z = 30;         %按照0.6s交织，则这改成30
end
kb = 6;
mb = 12;
nb = 18;
mod_odr = log2(M_ary);
rep_num = 2;

%瑞利衰落信道模型
fd = 1;
delay = [0 0.00209];   % 多径时延
pd = [0 0];  
ch = rayleighchan(1/Fs, fd, delay, pd);     % 生成信道模型
ch.dopplerspectrum = doppler.gaussian;
ch.storepathgains = 1;
if (fd~=0)
    ch.storehistory = 1;
end
ch.resetbeforefiltering = 1;    % 跳频需要每次复位信道
order = 4;
mp = 2;


bit_len = kb * z;
code_len = nb * z;
syb_len = rep_num * code_len / mod_odr;
tblen = 50;
rz_fg = 1;
sap_len = syb_len * sps;
sta_len = bit_len;

t_seq = 0:Ts:(sap_len+sps)*Ts-Ts;
tdt_tilt = t_seq/Tb;
ca_tilt = exp(1i*pi*h_mod*(M_ary-1)*tdt_tilt);

TREL = msk_trellis_gen(g_t, Ts, sps, h_mod, ca_tilt, fig_on);

snr_start = -4.6;
snr_step = 0.2;
snr_end = 20;

 
for snr = (snr_start:snr_step:snr_end)
    
    % 求解噪声功率
    eb_n0 = 10 ^ (snr/10.0);
    n0 = 1/((eb_n0)+1);
    sigma = n0 /2;
    
    be = 0;       %总的错误比特数
    fe = 0;       %总的错误帧数
    de_be = 0;    %累加可检测错误比特数
    un_be = 0;    %累加不可检测错误比特数
    de_fe = 0;    %累加不可检测错误帧数
    un_fe = 0;
    n_iter = 0;
    
    iter_flag=1;
    iter_num=0;
    
    frm_num = 0;%数据帧数
    gogo = 1;
    while gogo == 1
        bit_tx = randi([0 1], 1, bit_len);
        
        cod_tx = step(hl_enc, bit_tx.').';
        
        if (syb_len==5)
            cod_tx = [0 1 0 0 1];
        end
        
        
        %交织
        switch inter_tp
            case 0
                cod_tx_inter = cod_tx;
            case 1
                cod_tx_inter = matintrlv (cod_tx,inter_line,code_len/inter_line); %对信号序列的交织
            case 2
                cod_tx_inter = inter_prm(cod_tx, inter_p);
        end

        % 重复
        if (rep_num==2)
            cod_tx_rep = [cod_tx_inter, cod_tx_inter];
        else
            cod_tx_rep = od_tx_inter;
        end

        syb_tx_ori = 1 - 2*cod_tx_rep;
        
        pos_n = length(find(syb_tx_ori==1));
        if (mod(pos_n,2)==0)
            zf = -1;
            syb_tx = [syb_tx_ori, -1];
        else
            zf = 1;
            syb_tx = [syb_tx_ori, 1];
        end

        [dat_tx, phi_tx] = msk_mod(syb_tx, g_t, Ts, sps, h_mod);
        sig_pow = var(dat_tx);
        if (fig_on)
            freq_veiw( dat_tx, 1/Ts );
        end
        
        [ dat_rx, dat_rx_e, h_0 ] = chan_rf( dat_tx, ch, fd, snr, chann_flag, order, mp, sps, sig_pow, BW_bb);
        noise_bl = dat_rx - dat_rx_e;
        snr_real = 10*log10(sig_pow/var(noise_bl));
        
%         if (no_noise)
%             dat_rx = dat_tx;
%         else
%             dat_rx = awgn(dat_tx, snr-10*log10(sps/mod_odr)+0*10*log10(nb/kb), 'measured');
% %             dat_rx = awgn(dat_tx, snr-10*log10(sps/mod_odr)+0*10*log10(nb/kb));
%         end

        dat_rx_filt = conv(dat_rx, h_lpf);
        dat_rx_filt = dat_rx_filt(odr_lpf+1:odr_lpf+length(dat_rx));
%         figure; plot(real(dat_rx_filt)); hold on; plot(real(dat_rx_e), 'r');
        
        dat_rx_tilt = dat_rx_filt .* ca_tilt;
        
        llr_d = zeros(1, length(syb_tx));
        % 此处开始迭代     
        err_l = zeros(1, iter_max+1);
        while (iter_num<=iter_max && iter_flag~=0)

        
%             syb_rx = msk_dem_viterbi(dat_rx_tilt, TREL, tblen, rz_fg);
            [llr_rx, syb_rx] = msk_dem_logmap(dat_rx_tilt, TREL, llr_d, rz_fg);
            
            % 重复合并
            if (rep_num==2)
                llr_rx_dere = llr_rx(1:code_len) + llr_rx(code_len+1:2*code_len);
            else
                llr_rx_dere = llr_rx;
            end
            
            
            llr = (1/sqrt(2)) * llr_rx_dere(1:code_len) / (8*n0);
        
%             be_l = length(find(cod_tx~=syb_rx(1:syb_len)));
%             be = be + be_l;
%             if (be_l~=0)
%                 fe = fe + 1;
%             end
%             sta_len = syb_len;
        
            % 解交织
            switch inter_tp
                case 0
                    llr_deint = llr;
                case 1
                    llr_deint = matdeintrlv (llr,inter_line,code_len/inter_line);
                case 2
                    llr_deint = deint_prm(llr, inter_p);
            end

        
            % 译码
            [ de_be, un_be, de_fe, un_fe, n_iter, cod_rx, llr_dp, c_flag, be_now] = decode_mat(...
                hl_dec, llr_deint, bit_tx, mb, nb, z, de_be, un_be, de_fe, un_fe, n_iter, 0);
%             be = de_be + un_be;
%             fe = de_fe + un_fe;
            bit_rx = cod_rx(1:bit_len);
            bit_rx_pc = (1-sign(syb_rx)) / 2;
            llr_dp = llr_dp - llr_deint;
%             llr_dp = llr_dp*0.8^iter_num;
            % 判断译码正确否，如果不正确，对llr_dp做交织，交织后重新做解调
        
            err_l(iter_num+1) = length(find(bit_rx~=bit_tx));

%             be = be + err_l;
            if (err_l(iter_num+1)~=0)
                iter_flag=1;
                % 再次交织
                switch inter_tp
                    case 0
                        llr_dd = llr_dp;
                    case 1
                        llr_dd = matintrlv (llr_dp,inter_line,code_len/inter_line);
                    case 2
                        llr_dd = inter_prm(llr_dp, inter_p);
                end

%                 pos_nn = length(find(cod_rx==0));
%                 if (mod(pos_nn,2)==0)
%                     llr_d = [llr_dp, -100];
%                 else
%                     llr_d = [llr_dp, 100];
%                 end
    
                % 再次重复
                if (rep_num==2)
                    llr_d = [llr_dd, llr_dd, zf*100];
                else
                    llr_d = [llr_dd, zf*100];
                end
                
            else
                iter_flag=0;
            end
%             iter_flag=0;
%             iter_num=0;

            iter_num = iter_num + 1;

        end
        
        be = be + err_l(iter_num);
        if(err_l(iter_num)~=0)
            fe=fe+1;
        end
        iter_num=0;
        iter_flag=1;
        frm_num = frm_num + 1;
          
        if (mod(frm_num,10)==0)
            bit_total = sta_len * frm_num;
            ber = be / bit_total;
            fer = fe / frm_num;
            fprintf('snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
        end
    
        if (fe>=100 && frm_num>=400)
            gogo=0;
        end
    end
    bit_total = sta_len * frm_num;
    ber = be / bit_total;
    fer = fe / frm_num;
    fprintf('snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
    fprintf(fid, 'snr=%.1f frm_num=%d ber=%.8f fer=%.8f be=%d fe=%d\n', snr, frm_num, ber, fer, be, fe);
end

total_time = toc; % STOPWATCH OFF

